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Abstract 

Consider  the  vector  initial  value  problem  ey  = f(t,y,c),  y(0)  = y°(e) 
with  f (t,y ,0)  =*  F(t)y  + G(t)  for  a singular  matrix  F(t)  of  constant  rank 
with  stable  eigenvalues  and  zero  eigenvalues  having  simple  elementary  divis- 
ors. This  paper  shows  how  to  determine  the  unique  limiting  solution  when 
the  reduced  problem  FYQ  + G = 0 is  solvable  and  how  to  obtain  a uniform 
asymptotic  expansion  for  the  solution  on  finite  t intervals. 


This  work  was  supported  in  part  by  the  Office  of  Naval  Research  under  Con- 
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1 . Introduction 

We  wish  to  consider  the  vector  initial  value  problem 
O)  ey  = f(t,y,e),  y(0)  = y°(e) 

in  the  nearly  linear  situation  that 

(2)  f(t,y,0)  = F(t)y  + G(t). 

When  F(t)  is  stable  on  some  interval  0 £ t £ T,  it  is  a classical  result 
of  Tihonov  [17]  that  the  limiting  solution  of  the  initial  value  problem  (1) 
will  converge  to  the  solution  Y^  = -F  of  the  reduced  equation 

(3)  0 «»  F(t)YQ  + G(t) 

for  0 < t T as  the  small  positive  parameter  e tends  to  zero.  (The 
result  continues  to  be  valid  for  all  t > 0 if  YQ  decays  exponentially  to 
zero  as  t -*■  « (cf.  Hoppensteadt  [7]).)  The  solution  generally  converges 
nonuniformly  at  t » 0 since  we  cannot  expect  that  y^(0)  = Y^(0) . Assuming 
infinite  differentiability  of  f in  t and  y and  asymptotic  series  expan- 
sions in  e for  both  f and  y^,  we  can  Indeed  show  that  the  asymptotic 
solution  y(t,e)  of  (1)  is  of  the  form 

(4)  y(t,e)  = Y(t,e)  + JI(t  ,e) 

throughout  0 t £ T where  Y and  II  have  asymptotic  expansions  in  t 
and  the  terms  of  n all  tend  to  zero  as  the  stretched  (or  boundary  layer) 
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variable 

(5) 


t = t/e 


tends  to  infinity.  Thus,  the  limiting  asymptotic  solution  for  t > 0 is 
provided  by  the  outer  solution  Y(t,e),  and  the  nonuniform  convergence  of 
the  solution  from  y^(e)  to  Y(0,c)  near  t = 0 is  accomplished  through 
the  boundary  layer  correction  II(T,e). 

If  F(t)  has  unstable  eigenvalues,  we  cannot  expect  the  initial  value 
problem  to  have  bounded  solutions  as  e -*•  0 unless  we  restrict  the  initial 
values  to  a lower  dimensional  manifold  (cf.  Levin  [11]).  It  is  still  reason- 
able, however,  to  allow  arbitrary  initial  values  if  the  matrix  F(t)  is 
singular.  (The  usual  expansion  procedure  will  then  break  down.)  For  an 
M-vector  y,  let  us  assume 


A 


(HI)  the  M x M matrix  F(t)  has  <a  constant  rank  k,  0 <_  k < H,  for  all 
t in  0 < t <_  T,  its  nonzero  eigenvalues  all  have  negative  real  parts 
there,  and  its  null  space  is  spanned  by  M - k linearly  independent 
eigenvectors. 

Under  (HI),  we  shall  find  that  the  asymptotic  solution  of  (1)  - (2)  re- 
mains in  the  form  (4)  whenever  the  limiting  equation  (3)  is  consistent,  i.e. 

G CR(F),  the  range  of  F.  Because  F is  singular,  (3)  no  longer  has  a 
unique  solution.  Its  solution  is  determined  up  to  an  arbitrary  element  of 
N(P),  the  null  space  of  F,  so  additional  analysis  is  required  to  determine 
the  unique  limiting  solution  for  t > 0.  This  is  an  instance  of  a "singular 
singular-perturbation  problem"  where  the  reduced  problem  (obtained  when 


_ 
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e ■ 0)  has  an  infinity  of  solutions  in  regions  of  uniform  convergence. 

More  generally,  one  could  consider  singular  problems  (1)  where 
fy(t,y,0)  - 0 for  some  t values  along  some  solutions  y of  the  reduced 
equation  f(t,y,0)  = 0.  Difficult  turning  point  problems  would  then  be  faced 
(cf.  Wasow  [20]).  A somewhat  more  restricted  class  of  nonlinear  problems  is 
considered  in  O'Malley  and  Flaherty  [15].  Their  analysis  seeks  to  eventually 
develop  numerical  algorithms  for  such  problems.  Among  other  possible  appli- 
cations are  high  gain  feedback  control  systems  (cf.  O'Malley  [13]),  singular 
control  problems  (O'Malley  [14]),  and  certain  Cauchy  problems  in  Banach 
spaces  (Gordon  [6]) . We  note  that  the  structure  of  the  asymptotic  solutions 
will  generally  differ  considerably  if  the  elementary  divisors  of  F corres- 
ponding to  zero  eigenvalues  are  nonlinear.  An  example  is  provided  by 

I -1  0 0 . 

* F(t,e)y  for  F ■ l 0 0 -e  ) . Further,  if  F remains  non- 

'0-1  o' 

singular,  but  has  purely  Imaginary  eigenvalues,  rapidly  oscillating  solutions 
would  result  (cf.  Hoppensteadt  and  Miranker  [9]). 

2.  Preliminary  Linear  Algebra 

Under  hypothesis  (HI) , we  are  guaranteed  that  the  matrix  F can  be  put 
into  its  row-reduced  form  by  an  orthogonal  matrix  E(t),  i.e. 


where  U(t)  is  a k x M matrix  of  rank  k.  Indeed,  E is  guaranteed  by 
the  singular  value  decomposition  of  F,  i.e.  F - E'DH  for  orthogonal 
matrices  E and  H and  diagonal  D (cf.,  e.g.,  Stewart  [16]).  E can  be 
obtained  numerically  by  performing  a sequence  of  Householder  transformations 
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(cf.  Golub  [4]).  Moreover,  the  differentiability  of  E will  follow  that  of 
the  coefficients  under  the  constancy  of  rank  assumption  (cf.  Golub  and 
Pereyra  [5]). 

The  use  of  such  orthogonal  matrices  to  reduce  F to  echelon  form  seems 
more  convenient  than  the  diagonalization  procedures  which  have  been  more 
common  in  the  singular  perturbation  literature.  Among  other  differences,  it 
avoids  the  use  of  eigenvectors  of  F.  The  success  of  traditional  procedures 

on  closely  related  problems  is  demonstrated  in  Fife  [3],  Butuzov  and 
Vasil'eva  [1],  and  Vasil'eva  [19].  Alternative  analogous  methods  which  don't 
involve  orthogonal  matrices  may  be  necessary  for  nonlinear  problems  (cf. 
O’Malley  [14]). 

Writing 


where  is  k x M,  we  have  the  orthogonality  condition 


e2f  - 0. 


Furthermore,  the  orthogonality  of  E implies  that 


eie2  ‘ °-  Vi  ■ V Vi  ■ W ^ EiEi  + e;e2  ■ V 


These  last  relationships  imply  that 


P - E^E1  - P2  and  Q - E’E2  - Q2 


(10) 
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are  complementary  projections.  Indeed,  a direct  sum  decomposition  of  M-space 
results  with 


(11) 


R(Q)  - W(F’)  and  R(P)=R(F). 


shows  that  Q projects  into  W(F') 
are  both  M - k.  Finally,  (7)  and 

where  the  k x k matrix 
stable 

since  it  shares  the  k stable  eigenvalues  of  F.  Indeed,  the  purpose  of 
(HI)  is  to  provide  an  orthogonal  matrix  E such  that  EF  is  row-reduced 
with  EFE'  having  k stable  eigenvalues. 

3.  The  Transformation  Approach 

Making  the  one-one  transformation 


This  follows  since  F'Q  = (E2p)'E2  = 0 
and  the  dimensions  of  R(Q)  and  W(F') 

lhni  VE2\ 

(8)  show  that  EFE'  = I i 

\ 0 0 / 

(12)  S = Ej^FE^  is 


(13)  z ■ E(t)y, 
z will  satisfy 

(14)  ez'  - EFE'z  + EG  + e[EE'z  + i E(f(t,E'z,E)  - fU.E'z.O))] . 


Setting 

(15) 
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then,  provides  a transformed  singular  perturbation  problem 


(16) 


(ez|  = + E1FE12z2  + EjG  + eg1(z1,z2,t,e) 

EZ2  = E2G  + e82(zl’Z2,t,E^ 


1 

where  g2  = E2E'z  + — E2(f (t,E'z,e)  - f(t,E'z,0)).  The  asymptotic  behavior 
of  solutions  to  (16)  differs  considerably  in  the  two  cases  E2G  = 0 and 
E2G  + 0.  Since  EjG  =0  if  and  only  if  QG  = 0,  Q projecting  onto  W(F'), 
we're  guaranteed  that  EjG  =0  if  and  only  if  the  reduced  problem  (3)  is 
consistent.  Taking 


(17) 


E2G  = 0 throughout  [0,T], 


then,  reduces  (16)  to  a "nonsingular  singular-perturbation  problem."  For  it, 
classical  singular  perturbation  theory  (cf.,  e.g.,  Hoppensteadt  [8]  or 
O'Malley  fl2j)  implies  that  (16)  - (17)  has  a limiting  solution  for  t > 0 


which  is  the  unique  solution 


'10 


'20 


of  the  reduced  problem 


0 - SZ10  + EjPEJZjd  + ElG 


(18) 


^20  " ®2^Z10,Z20’*' m Eo(0)y  (0)» 


20' 


presuming  that  the  corresponding  nonlinear  initial  value  problem 


(19) 


'20 


g2(-S"1E1(FE^Z20  + G),Z20,t,0), 


z2o(0)  ■ E2(0)y  (0) 
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has  a solution  for  0 <_  t <_  T.  Furthermore,  that  theory  also  shows  that  the 
nonuniform  convergence  of  z near  t = 0 is  determined  as  the  decaying 
solution  of  the  boundary  layer  problem 


{ dT 


10 


(20) 


dr 

dT 


S(0)T 


10’ 


Tl0(°) 


E1(0)yu(0)  - Z1(0) 


20 


\ dr 


Thus,  T^q(t)  = e ^ Tiq(0)  and  T2q(0)  = 0 for  t ^ 0.  The  bifurcation 
character  (cf.,  e.g.,  Cesari  [2])  of  these  limiting  problems  makes  them  dif- 
fer from  the  classical  ones.  Here,  the  reduced  problem  (18)  consists  of  an 
algebraic  problem  for  a k-vector  and  an  initial  value  problem  for  an  M - k 
vector,  while  the  boundary  layer  problem  (20)  involves  a k dimensional 
initial  value  problem  and  an  M - k dimensional  algebraic  problem. 


4.  A Direct  Solution 

Since  the  boundary  layer  correction  n(r,e)  will  be  asymptotically 
negligible  for  each  t > 0,  the  outer  solution 


(21)  Y(t,e)  -v  l Y (t)eJ 

j-0  3 

must  satisfy  the  original  nearly  linear  system  (1)  as  a power  series  in  e. 
Rewriting  the  equation  as 

(22)  F(t)Y  - -G(t)  - [f (t,Y,e)  - f(t,Y,0)]  + cY 


and  substituting  in  the  expansion  (21)  Implies  that  the  terms  of  (21)  must 


successively  satisfy  the  algebraic  equations 


(23)  F(t)Y  = Cj_1(t) 

where 


(24) 


/ and 


-G, 


fe(t,Y0,0), 


V 


fye(t’V0)Yk  + nk-l(t)’ 


k > 1 


for  a known  function  ^ of  Y^,  Y^,  ...»  Y^  Assuming  consistency  of 
(3),  G c R(F)  so,  by  (8),  QG  = 0.  Consistency  of  (23),  by  the  Fredholm 
alternative,  requires  ^ = 0,  j >_  1,  so  we  obtain  the  singular  differ- 

ential equations 


(25) 


” Qy0  - Q*«Ct,Y0,o) 

and 

Q\  - Qfye(t,Y0,0)Vk  - k > 1. 


Taking  (23)  and  (25)  together  nearly  allows  us  to  obtain  the  outer  solu- 
tion Y termwise.  Specifically,  multiplying  (23)  by  implies  that 

E^F(E^E^  + Q)Yj  ■ E^Cj  j ^ 0,  so  the  invertibility  of  S * E^FE^  and  the 
definition  of  P imply  that 


(26) 


PYj  = -AF(QYj)  + ACj_1 


where  A - E^S  ^E^.  Thus  PY^  and  thereby  are  completely  determined  as 
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linear  functions  of  QY^ . Putting  (25)  and  (26)  together,  then,  implies  the 
differential  equations 

(QY0)'  = Q(B(QYq)  - AG)  + Qf£(t,8(QY0)  - AG.O) 

(27)  < and 

k (QYk)'  = [Q  + Qfye(t,Y0,0)][8(QYk)  + A5kl)  - Qnk_1 , k > 1 
since  (26)  implies  that 

(28)  Yj  = 8(QY^)  + A? 

for  the  projection  8 = I - Af.  Note  that  the  representation  (4)  implies  the 
yet  unspecified  initial  values 

(29)  Q(0)Yk(0)  = Q(0)yk  - Q(0)IIk(0) , k > 0 

needed  to  completely  determine  the  outer  solution  Y(t,e). 

If  we  anticipate  that  Q(0)nQ(0)  = 0 (paralleling  the  result  that 
T2q(0)  “0  in  (20)),  it  will  follow  that  the  limiting  outer  solution  is  given 
by 

(30)  YQ  - 8(QYQ)  - AG 

where  QYq.  the  projection  of  the  limiting  solution  onto  M(F'),  satisfies 
the  nonlinear  initial  value  problem 


(QYq)  = Q(8(QY0)  - AG)  + Qfe(t,B(QY0)  - AG,0) 


Q(0)Yq(0)  = Q(0)y  (0). 


We  shall  assume  that  the  solution  to  (31)  exists  throughout  0 <_  t <_  T (in 
analogy  to  our  existence  assumption  for  (19)).  We  note  that  (31)  shows  that 
the  limiting  solution  of  (1)  - (2)  is  determined  by  a dynamical  system  in  an 
M - k dimensional  space.  Further,  we  observe  that  unique  solutions  QY^  and 
Y^  for  k > 0 follow  from  the  linear  differential  equations  of  (27)  and  the 
algebraic  relation  (28). 


By  the  form  of  (4),  the  boundary  layer  correction  IT  must  satisfy  the 


system 


dll  dy  dY 

= f _ c = 


dx  dt 


e — = f (ex ,Y (ex , e)  + n(x,e),e)  - f (ex,Y(ex,e) ,e) 


and  the  terms  of  its  expansion 


II(x,e)  n.  T.  II.  (x)e^ 
j=0  J 


should  tend  to  zero  as  x tends  to  infinity.  In  the  quasilinear  situation 
when  (2)  holds,  equating  coefficients  successively  in  (32)  implies  the  linear 
system 


' r(°>"i  + Vi't) 


where  ^ is  a linear  combination  of  the  preceding  coefficients  JT^,  t < j 
with  polynomial  coefficients  in  x,  and  c , = 0.  Since  QF  = 0,  it  follows 
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that  — (Q(0)I1,(t))  = Q(0)C.  , (t).  Thus,  the  decay  requirement  implies  that 

dT  j j-1 


Q(0)nj(r)  = -Q(0)  | cj_1(s)ds. 


We  note,  in  particular,  that  (35)  allows  the  termwise  determination  of  the 
initial  values  (29)  needed  to  fix  the  outer  solution.  It  remains  to  succes- 
sively find  the  P(0)IK(t)'s.  Multiplying  (34)  by  E^(0)  implies  that 
^ (E1(0)nj)  = S(0)(E1(0)nj)  + E1(0)F(0)(Q(0)nj)  + E1(0)cj_1.  Using  (35)  and 

Iu  = E'(0)E,(0)  + Q(0)  implies  the  decaying  solutions 
Mil 


(36)  jyx)  = E|(0)eS(0)TE1(0)n.(0)  - Q(0) 


(s)ds 


,T  f°° 

!’(0)  eS(0)(T_r)E1(0)[Cjl(r)  - F(0)Q(0)  J ^ _1(s)ds]dr. 


Thus,  the  boundary  layer  terms  are  determined  successively  as  exponentially 
decaying  vectors  up  to  the  initial  value  E^(0)II(0,e) , just  as  the  outer  ex- 
pansion coefficients  are  determined  up  to  Q(0)Y(0,e).  Since  the  limiting 
outer  solution  is  known  from  (30)  - (31)  and  higher  order  terms  follow  succes- 
sively from  (27)  - (29),  the  unique  higher  order  boundary  layer  correction 
terms  follow  from  (36)  with 


E1(0)nj(0)  - E1(0)(yj  - Yj (0)) , j > 0. 


For  j ” 0,  we'll  have 


(38)  nQ(T)  = E|(0)eS(0)TE1(0)(y°(0)  - 8(0)Q(0)y°(0)  - A(0)G(0)). 


Since  QEj  ■ 0,  it  follows  that  11^  lies  within  the  k dimensional 


I 
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R(F)  = R(P).  The  function  of  this  zero  order  boundary  layer  correction  term, 
then,  is  to  "instantly"  transfer  from  the  prescribed  initial  vector  y^(0) 
to  Yq(0+)  where  the  reduced  equation  (3)  is  satisfied.  Finally,  we  observe 
that  the  expansion  (A)  which  we've  formally  generated  will,  of  course,  agree 
with  that  obtained  by  our  earlier  transformation  procedure  and  is  justified 
by  it. 

Summarizing,  we  have 
Theorem 

Under  hypothesis  (HI);  (H2),  that  (3)  is  consistent;  and  (H3) , that  the 
solution  (QYq)  £f  (3A)  exists  for  0 _<  t _<  T,  the  unique  solution  of  the 
initial  value  problem  (1)  - (2)  i£  asymptotically  of  the  form  (A)  uniformly 
in  Oj^t^T  as  e -►  0. 

5.  The  Inconsistent  Problem 

If  E2G  i 0,  the  transformed  system  (16)  cannot  be  easily  solved  as 
(16)  - (17)  was.  Then,  the  limiting  system  (3)  has  no  solutions  and  we  cannot 
expect  to  have  an  outer  solution  of  the  form  (21). 

By  variation  of  parameters,  we  know  that  solution  of  (1)  - (2)  must  sat- 
isfy the  integral  equation 

(39)  y(t,c)  - ^(t ,e)y°(e)  + 7 f y'(t,e)y'"1(s,e){G(t)  + f(t,y,e)  - f(t,y,0)}dt 

E J0 

where  V is  the  fundamental  matrix  solution  of  the  linear  system 

(40)  y « i F(t)y,  y(0,e)  - 1^ 

Such  fundamental  matrices  can  be  constructed  following  the  work  of  Turrittin 
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and  earlier  authors  (cf.  Turrlttin  [18]).  Under  natural  assumptions,  we  can 
show  that  ||y(t,e)y  ^(s,e)||  K for  0 <_  s,t  T for  some  constant  K 
provided  no  eigenvalues  of  F are  unstable.  When  F is  stable,  we  instead 
have  a decaying  bound  Ke  S^E  and  this  allows  the  corresponding  outer 
solution  for  (1)  - (2)  to  be  bounded.  When  F has  eigenvalues  with  zero  real 
parts,  other  estimates  are  appropriate  (cf.  Hoppensteadt  and  Miranker  [9]). 

For 


(41)  f(t,y,e)  = F(t,e)y  + G(t,e) 

linear  in  y,  (39)  makes  it  natural  to  seek  the  solution  to  (1)  in  the  form 

(42)  y(t,c)  - ^ Z(t,e)  + n(x,e) 

CD 

when  hypothesis  (HI)  holds  for  F with  (F(t,e) ,C(t,e))  E (F  (t),G  (t))e^. 

j-0  J j 

Here  the  outer  solution  Z(t,c)  and  the  boundary  layer  correction  II  have 
power  series  expansions  in  c and  n -*■  0 as  t -*•  «>.  Now  the  limiting  solu- 
tion ~ Zg(t)  for  t > 0 satisfies  the  homogeneous  system  F^(t)Z^  - 0. 
Applying  the  Fredholm  alternative  to  the  equation  for  F^Z^  requires  that 
Q(Zq  - FjZq  - Gq)  ■ 0 where  Q is  the  projection  onto  N(F^).  Using  our 
earlier  notation,  then,  we  have  the  unbounded  limiting  solution  Z^(t)  with 

(43)  ZQ  - 8(QZq) 
where  QZq  satisfies  the  linear  system 

(44)  (QZq) * - (Q  + QF1)8(QZ0)  + QGq,  Q(0)Zq(0)  - 0. 


_ 
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The  boundary  layer  correction  is  calculated  as  before.  When  QG^  = 0, 

Zg  = 0 and  the  solution  reduces  to  the  form  (4)  obtained  for  the  consistent 
problem. 

If  f is  nonlinear  in  y,  the  expansion  (42)  cannot  generally  be  used, 
k lk 

(A  term  like  y would  provide  a term  — in  the  outer  region.)  For 

e 

polynomial  dependence  on  y,  an  appropriate  scheme  might  be  devised  for  the 
inconsistent  problem  though  we  shall  not  do  so  (cf.  Kersten  [10],  however, 
for  a related  discussion  concerning  special  boundary  value  problems). 
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